#!/usr/bin/env julia
# Capture summary statistics for a radiance BUFR file using the native ecCodes
# reader and store them in JSON for regression comparisons.
#
# Usage:
#     julia --project=. scripts/julia/capture_radiance_stats.jl <bufr_file> <output_json>
#
# The script infers instrument metadata from the decoded observation structure and
# records channel means, standard deviations, min/max counts, and observation
# geometry ranges. The emitted JSON is suitable for baseline comparisons against
# Fortran GSI diagnostics or prior Julia runs.
using JSON3
using Dates
using Statistics
using GSICoreAnalysis
using GSICoreAnalysis.DataIO
using GSICoreAnalysis.DataIO.RadianceBUFRReader
import GSICoreAnalysis.DataIO.RadianceBUFRReader: summarize_radiance_observations

function summarize(bufr_file::String; max_obs::Union{Int,Nothing}=nothing)
    config = RadianceBUFRConfig(verbose=false, max_obs=max_obs)
    obs = read_radiance_bufr(bufr_file; config=config)

    summary = summarize_radiance_observations(obs)

    nchans, nobs = size(obs.brightness_temp)
    channel_stats = Dict{String,Any}()
    for (idx, ch) in enumerate(obs.channels)
        bt = obs.brightness_temp[idx, :]
        valid = filter(!isnan, bt)
        channel_stats[string(ch)] = Dict(
            "mean" => isempty(valid) ? missing : mean(valid),
            "std" => isempty(valid) ? missing : std(valid),
            "valid_count" => length(valid),
        )
    end

    geom = Dict(
        "latitude_min" => minimum(obs.latitude),
        "latitude_max" => maximum(obs.latitude),
        "longitude_min" => minimum(obs.longitude),
        "longitude_max" => maximum(obs.longitude),
        "zenith_min" => minimum(obs.zenith_angle),
        "zenith_max" => maximum(obs.zenith_angle),
    )

    return Dict(
        "source_file" => abspath(bufr_file),
        "generated_at" => string(Dates.now()),
        "instrument" => obs.instrument,
        "satellite" => obs.satellite_id,
        "n_channels" => nchans,
        "n_observations" => nobs,
        "channel_statistics" => channel_stats,
        "geometry" => geom,
        "metadata" => summary,
    )
end

function main()
    if length(ARGS) < 2 || length(ARGS) > 3
        println(stderr, "Usage: julia --project=. scripts/julia/capture_radiance_stats.jl <bufr_file> <output_json> [max_obs]")
        exit(1)
    end

    bufr_file = ARGS[1]
    output_json = ARGS[2]
    max_obs = if length(ARGS) == 3
        try
            parse(Int, ARGS[3])
        catch
            println(stderr, "ERROR: max_obs must be an integer")
            exit(3)
        end
    else
        nothing
    end
    if !isfile(bufr_file)
        println(stderr, "ERROR: BUFR file not found at $bufr_file")
        exit(2)
    end

    stats = summarize(bufr_file; max_obs=max_obs)
    open(output_json, "w") do io
        JSON3.write(io, stats; indent=2, allow_inf=true)
    end
    println("Wrote radiance statistics to $output_json")
end

main()
